#!/usr/bin/perl

#Written By: Matthieu Chartier
#Adapted By: Aaron Steele

#Build translation table hash
open TRANSL_TBL, "<files/transl_tbl" or die "Cant open transl_tbl file";
foreach(<TRANSL_TBL>){
    @elements=split(/,/, $_);
}
foreach(@elements){
    $_=~s/^\s+//;
    $transl_tbl{$1}=$2 if($_=~/'([a-z0-9\*]+)'=\>'([a-z0-9\*]+)'/i);
}
print "\nDone building translation table\n";

#Open codon usage file and store the table as an array
open CUT, "files/cut" or die "\nCan't open file codon usage table file\n";
foreach(<CUT>){
    push @cut, $_;
}
close CUT;

unless($start_at_s3 eq "y"){
    print "\n--------------- STAGE 1 ---------------\n";
    print "Retrieving info for each Pfam\n";
    @pfam_auto_files_list = glob "pfam_files/*";
    foreach(@pfam_auto_files_list){
	chomp($pfam_file = $_);
	$pfam_file =~ s/pfam_files\///;
	print "\n$pfam_file\n";
	
	open CONTENT, "pfam_files/$pfam_file" or die "\nCan't open PFAM alignment file $pfam_file\n";
	foreach(<CONTENT>){
	    push @pfam_file, $_;
	}
	$nb_index = $#pfam_file;
	
	#Open file in which to write data
	open ALIGN, ">stage1/pfam_$pfam_file" or die "Cant open output file in pfam_output";
	
	#Search the pfam alignment for the line at which alignment starts
	for($i=0; $i<=$nb_index; $i++){
	    if($pfam_file[$i] =~ /^[^#]/){
		$index = $i;
		last;
	    }
	}
	
	###########################################################
	#-----| FOR EACH PROTEIN IN THE PFAM ALIGNMENT FILE |-----#
	#-----|-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-|-----#
	###########################################################
	
	$protein_count=0;
	$protein_stored=0;
	$protein_no_cu=0;
	$protein_no_ebi=0;
	$protein_aa_not_matched=0;
	$protein_wrong_codon_nb=0;
	$protein_nuc_no_match=0;
	
      PROTEIN: for($z=$index; $z<=$nb_index; $z++){
	  if($pfam_file[$z] =~ /^([A-Z0-9]+_?[0-9A-Za-z]*)\//){
	      $uniprot_id = $1;
	      $protein_count++;
	      
	      #Get content from the uniprot website for this protein
	      $uniprot_url = "http://www.uniprot.org/uniprot/$uniprot_id";
	      use LWP::Simple;
	      @uniprot_content = ();
	      @split_uniprot_content = ();
	      @uniprot_content = get $uniprot_url or die "\nCan't get content from $uniprot_url";	
	      @split_uniprot_content = split(/\n/, "@uniprot_content");
	      
	      #Get organism's taxonomic ID	
	      foreach(@split_uniprot_content) {
		  if($_ =~ /http:\/\/www.ncbi.nlm.nih.gov\/Taxonomy\/Browser\/wwwtax.cgi\?lvl=0&amp;id=(\d+)"/i) {
		      $org_id = $1;
		      last;
		  }
	      }
	      
	      #Get Uniprot identifier
	      foreach(@split_uniprot_content){
		  $uniprot_id = $1 if($_ =~ /\<strong\>\s+([0-9a-z]+)\<\/strong\>\s+\($uniprot_id\)/i);
	      }
	      
	      $matched="0";
	      #Check if codon usage table of this organism is available
	      foreach(@cut){
		  if(/=($org_id)\|/){
		      $matched="1";
		      #print "\nCodon usage found";
		      last;
		  }
	      }
	      if($matched=="0"){ #Error #1
		  #print "\nCodon usage not available\n";
		  $protein_no_cu++;
		  push @protein_no_cu, $uniprot_id;
		  next PROTEIN;
	      }
	      
	      if($matched == "1"){
		  #Get pfam_ori and pfam_aa
		  chomp($pfam_align=$pfam_file[$z]);
		  $position=rindex($pfam_align, " ");
		  $pfam_aa_ori=substr $pfam_align, ($position + 1);
		  $pfam_aa_ori=~tr/a-z/A-Z/;
		  $pfam_aa=substr $pfam_align, ($position + 1);
		  $pfam_aa=~s/\.//g;
		  $pfam_aa=~s/-//g;
		  $pfam_aa=~tr/a-z/A-Z/;
		  
		  #Find ebi url from uniprot website
		  $ebi_url="";
		  foreach(@split_uniprot_content){
		      if($_=~/\<a\s+class="embl_cds"\s+href="(http:\/\/www\.ebi\.ac\.uk\/ena\/data\/view\/[\w]+)"\s+/i){
			  $ebi_url=$1;
			  last;
		      }
		  }
		  
		  if($ebi_url eq ""){ #Error #2
		      #print "\nCouldn't find cds from ebi.";
		      $protein_no_ebi++;
		      push @protein_no_ebi, $uniprot_id;
		      next PROTEIN;
		  }
		  $ebi_url.="&display=text";
		  #print "\n$ebi_url";
		  
		  ############################################################
		  #########-------Fetch nucseq from EBI website-------########
		  
		  use LWP::Simple;
		  @content = get("$ebi_url") or die "\nCan't get content from $ebi_url\n";
		  @split_content = split(/\n/, "@content");
		  $i_max = $#split_content;
		  
		  ##########################################
		  #-------GET TRANSLATION TBL NUMBER-------#
		  $transl_tbl="";
		  foreach(@split_content){
		      if(/\/transl_table=([\d]+)/){
			  $transl_tbl=$1;
		      }
		  }
		  $transl_tbl="1" if($transl_tbl eq ""); #If transl tbl nb is not mentionned, its because transl tbl 1 is used
		  #print "\nTransl_tbl: $transl_tbl";
		  
		  #########################################
		  #------------GET TRANSLATION------------#
		  
		  #Search line at which translation starts
		  for($i=0; $i<=$i_max; $i++){
		      if($split_content[$i] =~ /\/translation="/i) {
			  $translation_line_start = $i;
			  last;
		      }
		  }
		  
		  #First, check if translation starts and stops at the same line
		  $solved = "2";
		  for($i=$translation_line_start; $i<=$i_max; $i++){
		      if($split_content[$i] =~ /\/translation="([\w]+)"/i){
			  $translation = $1;
			  #print ALIGN "\nncbi_aa=$translation"; #Write translation into the file
			  $solved = "1";
			  last;
		      }
		  }
		  
		  #If not, find where translation stops
		  if($solved != "1"){
		      for($i=($translation_line_start + 1); $i<$i_max; $i++){
			  if($split_content[$i] =~ /"/){
			      $translation_line_stop = ($i);
			      last;
			  }
		      }
		      
		      #Store the first line into $translation
		      $position = (index($split_content[$translation_line_start], "\"")) + 1;
		      $firstline = substr($split_content[$translation_line_start], $position);
		      $translation = $firstline;		
		      
		      #Add the middle lines to $translation
		      $i = (($translation_line_start) + 1);
		      while($i < $translation_line_stop){
			  $line = substr($split_content[$i], 2);
			  if($line =~ /\w/i){
			      $translation .= "$&$'";
			  }
			  $i++
		      }
		      
		      #Add the last line to $translation
		      $lastline = substr($split_content[$translation_line_stop], 2);
		      if($lastline =~ /\s([\w]+)"/i){
			  $translation .= "$1";
			  #print "\nTranslation stored.";
		      }	       
		  }
		  
		  #see if the pfam_aa can be matched in the translation if yes get the dna sequence
		  if($translation =~ /($pfam_aa)/){
		      #print "\nFound pfam_aa in translation.";
		      
		      ##########################################
		      #------------GET DNA SEQUENCE------------#
		      
		      #Search line at which DNA sequence starts
		      for($i=0; $i<=$i_max; $i++){
			  if($split_content[$i] =~ /SQ\s+Sequence/i) {
			      $seq_line_start = ($i +1);
			      last;
			  }
		      }
		      
		      #Search line at which DNA sequence stops
		      for($i=$seq_line_start; $i<=$i_max; $i++){
			  if($split_content[$i] =~ /\/\//){
			      $seq_line_stop = ($i);
			      last;
			  }
		      }
		      
		      #Store DNA sequence into variable
		      $i = $seq_line_start;
		      while($i < $seq_line_stop){ #for each line, do the following:
			  $_ = $split_content[$i];
			  s/\s+//g; #remove spaces
			  s/\d+//g; #remove nucleotide position
			  $dna_seq .= $_; #write the line in the file
			  $i++
		      }
		      $dna_seq = uc($dna_seq);
		      
		      #Find the nucleotide sequence in dna_seq that codes for the pfam_aa
		      @coding_seq = ();
		      $nuc_tot = length($dna_seq);
		      $pfam_aa_length = length($pfam_aa);
		      $i = 0; #reset position to 0 in dna_seq
		      $pfam_aa_pos = 0; #reset position to 0 in pfam_aa seq
		      $found_match_pos = 0;
		      while($i < $nuc_tot){ #parse each codon of dna_seq until it finds the seq that matches the PFAM_AA
			  $codon = substr($dna_seq,$i,3); #store codon
			  $aa_from_nuc_seq = $transl_tbl{$codon.$transl_tbl};
			  $aa = substr($pfam_aa,$pfam_aa_pos,1); #store aa
			  if($aa_from_nuc_seq eq $aa){ #check if codon gives same aa as the aa from pfam_aa
			      $found_match_pos = $i if($found_match_pos == 0);
			      push @coding_seq, $codon;
			      $pfam_aa_pos++;
			      if($pfam_aa_pos == $pfam_aa_length){ #stop if we reached the end of the aa seq
				  $nuc_seq_matched=1;
				  last;
			      }
			      $i += 3; #continue on to next codon
			  }
			  else{ #if there's no match: go back to begining of aa seq and clear coding_seq array
			      if($found_match_pos !=0){
				  $i = $found_match_pos+3;
				  $found_match_pos = 0;
			      }
			      else{
				  $i += 3;
			      }
			      $pfam_aa_pos = 0;
			      @coding_seq = ();
			  }
		      }
		      
		      #Store nucleotide sequence in a string
		      foreach(@coding_seq){
			  $coding_seq_string .= $_; 
		      }
		      
		      #Check if we haven't found the nucleotide sequence that codes for the pfam_aa chain
		      if($nuc_seq_matched == 0){ #Error #4
			  #print "\nCouldn't find nucleotide sequence corresponding to aa seq!";
			  $protein_nuc_no_match++;
			  push @protein_nuc_no_match, $uniprot_id;
			  next PROTEIN;
		      }
		      
		      $nb_codons_seq = (length($coding_seq_string))/3;
		      
		      #If nb of codons = nb of aa in the pfam alignment store all info into the file
		      if(($pfam_aa_length == $nb_codons_seq) && $nuc_seq_matched ==1){
			  $protein_stored++;
			  
			  $coding_seq_string = uc($coding_seq_string);
			  #print "\nNb of codons matched nb of aa!\n";
			  
			  #Write data into DATA file
			  print ALIGN "u_id=$uniprot_id//";
			  print ALIGN "\norg_id=$org_id//";
			  print ALIGN "\npfam_ori=$pfam_aa_ori//";
			  print ALIGN "\npfam_aa=$pfam_aa//";
			  print ALIGN "\ntransl_tbl=$transl_tbl//";
			  print ALIGN "\ncoding_seq=$coding_seq_string//\n";
		      }
		      elsif(($pfam_aa_length != $nb_codons_seq) && $nuc_seq_matched ==1){ #Error #5
			  #print "\nNb of codons DIDN'T match nb of aa!";
			  $protein_wrong_codon_nb++;
			  push @protein_wrong_codon_nb, $uniprot_id;
			  next PROTEIN;
		      }
		  }
		  else{ #Error #3
		      #print "\nCounld't find pfam_aa in translation!";
		      $protein_aa_not_matched++;
		      push @protein_aa_not_matched, $uniprot_id;
		      next PROTEIN;
		  }
	      }
	  }
	  $nuc_seq_matched=0;
	  $matched = 0;
	  $translation = "";	  
	  $translation_line_start = "";
	  $translation_line_stop = "";
	  $dna_seq = "";
	  $seq_line_start = "";
	  $seq_line_stop = "";
	  $coding_seq_string = "";
	  $pfam_aa_length = "";
	  $ebi_url = "";
	  $nb_codons_seq = "";
	  $i_max = "";
      }
	close CONTENT;
	close ALIGN;
	
	###------DATA EXTRACTION STATISTICS------###
	#print "\n\n------------------------------------";
	#print "\n-----DATA EXTRACTION STATISTICS-----";
	#print "\n------------------------------------";
	print "$protein_stored proteins over $protein_count passed stage 1. See extraction log file for details.";
	
	if($protein_stored > 0){
	    open LOG_ERRORS, ">results/$store_res_dir/extraction_log" or die "Cant open error log file";
	    print LOG_ERRORS "Proteins with no codon usage info\n" if($#protein_no_cu>=0);
	    print LOG_ERRORS "$_ $pfam_file\n" foreach(@protein_no_cu);
	    
	    print LOG_ERRORS "\nProteins for which the nucleotide sequence couldn't be retrived\n" if($#protein_no_ebi>=0);
	    print LOG_ERRORS "$_ $pfam_file\n" foreach(@protein_no_ebi);
	    
	    print LOG_ERRORS "\nProteins for which the aa sequence in the pfam alignment couldn't be matched with the one on the uniprot website\n" if($#protein_aa_not_matched>=0);
	    print LOG_ERRORS "$_ $pfam_file\n" foreach(@protein_aa_not_matched);
	    
	    print LOG_ERRORS "\nProteins with wrong number of codons\n" if($#protein_wrong_codon_nb>=0);
	    print LOG_ERRORS "$_ $pfam_file\n" foreach(@protein_wrong_codon_nb);
	    
	    print LOG_ERRORS "\nProteins for which the aa seq from the alignement wasn't matched with the nucleotide sequence\n" if($#protein_nuc_no_match>=0);
	    print LOG_ERRORS "$_ $pfam_file\n" foreach(@protein_nuc_no_match);
	    
	    close LOG_ERRORS;
	}
	
	$file_path = "stage1/pfam_".$pfam_file;
	unlink($file_path) if($protein_stored == 0);
	
	@pfam_file=();
	@protein_no_cu=();
	@protein_no_ebi=();
	@protein_aa_not_matched=();
	@protein_wrong_codon_nb=();
	@protein_nuc_no_match=();
	$protein_stored=0;
	
    }
    print "\n\nAll Nucleotide Sequences acquired.......Exiting!\n";
    exit 0;
}	

